module big_int_mod
    implicit none
    integer,parameter::maxn=1000
    integer(kind=4),parameter::maxnum=9
    type big_int
        integer(kind=4)::n(maxn)=0
        integer(kind=4)::s=1
        integer::h
    contains
        generic::assignment(=)=>init_big_int_char,init_big_int_int
        generic::operator(+)=>add
        generic::operator(-)=>sub
        generic::operator(*)=>mul_big_big,mul_big_small
        generic::operator(/)=>div
        generic::operator(==)=>compare_eq
        generic::operator(>)=>compare_gt
        generic::operator(<)=>compare_lt
        generic::operator(>=)=>compare_ge
        generic::operator(<=)=>compare_le
        generic::operator(/=)=>compare_neq
        procedure,pass::init_big_int_char,init_big_int_int
        procedure,pass::write_big_int
        procedure,pass::add,sub,mul_big_big,mul_big_small,div
        procedure,pass::compare_eq,compare_gt,compare_lt,compare_neq,compare_le,compare_ge
    end type
    interface mod
        module procedure mod_big
    end interface mod
    type(big_int),parameter::big_int_zero=big_int(0,1,1)
contains


    function compare_gt_pos(this,b)result(flag)
        type(big_int),intent(in)::this,b
        integer::ha,hb
        logical::flag
        integer::i
        ha=this%h
        hb=b%h
        flag=compare_eq_pos(this, b)
        if(flag)then
            flag=.false.
            return
        end if
        if(ha>hb)then
            flag=.true.
        elseif(ha<hb)then
            flag=.false.
        else
            do i=ha,1,-1
                if(this%n(i)>b%n(i))then
                    flag=.true.
                    exit
                elseif(this%n(i)<b%n(i))then
                    flag=.false.
                    exit
                end if
            end do
        end if
    end function compare_gt_pos

    function compare_gt(this,b)result(flag)
        class(big_int),intent(in)::this,b
        logical::flag
        integer::sa,sb
        flag=(this==b)
        if(flag)then
            flag=.false.
            return
        end if
        sa=this%s
        sb=b%s
        if(sa>sb)then
            flag=.true.
        elseif(sa<sb)then
            flag=.false.
        else
            flag=compare_gt_pos(this, b)
            if(sa==-1)then
                flag=.not.flag
            end if
        end if
    end function compare_gt

    function compare_lt(this,b)result(flag)
        class(big_int),intent(in)::this,b
        logical::flag
        integer::sa,sb
        flag=(this==b)
        if(flag)then
            flag=.false.
            return
        end if
        sa=this%s
        sb=b%s
        if(sa>sb)then
            flag=.false.
        elseif(sa<sb)then
            flag=.true.
        else
            flag=compare_gt_pos(this, b)
            if(sa==1)then
                flag=.not.flag
            end if
        end if
    end function compare_lt

    function compare_eq_pos(this,b)result(flag)
        type(big_int),intent(in)::this,b
        integer::ha,hb
        logical::flag
        ha=this%h
        hb=b%h
        flag=.false.
        if(ha==hb)then
            if(all(this%n(1:ha)==b%n(1:hb)))flag=.true.
        end if
    end function compare_eq_pos

    function compare_eq(this,b)result(flag)
        class(big_int),intent(in)::this,b
        logical::flag
        if(this%s==b%s)then
            flag=compare_eq_pos(this, b)
        else
            flag=.false.
        end if
    end function compare_eq

    function compare_neq(this,b)result(flag)
        class(big_int),intent(in)::this,b
        logical::flag
        flag=.not. (this==b)
    end function compare_neq

    function compare_ge(this,b)result(flag)
        class(big_int),intent(in)::this,b
        logical::flag
        flag=.not. (this<b)
    end function compare_ge

    function compare_le(this,b)result(flag)
        class(big_int),intent(in)::this,b
        logical::flag
        flag=.not. (this>b)
    end function compare_le

    function add_pos(this,b)result(c)
        type(big_int),intent(in)::this,b
        type(big_int)::c
        integer::ha
        integer(kind=4)::a,carry
        integer::i
        ha=max(this%h,b%h)
        carry=0
        c%s=1
        do i=1,ha
            a=this%n(i)+b%n(i)+carry
            if(a>maxnum)then
                c%n(i)=a-maxnum-1
                carry=1
            else
                c%n(i)=a
                carry=0
            end if
        end do
        if(carry==1)then
            c%n(ha+1)=1
            c%h=ha+1
        else
            c%h=ha
        end if
    end function add_pos

    function sub_pos(this,b)result(c)
        type(big_int),intent(in)::this,b
        type(big_int)::c
        integer::ha,hb
        integer(kind=4)::a,carry
        logical::flag
        integer::i
        ha=this%h
        hb=b%h
        flag=.true.
        if(ha>hb)then
            flag=.true.
        elseif(ha<hb)then
            flag=.false.
        else
            do i=ha,1,-1
                if(this%n(i)>b%n(i))then
                    flag=.true.
                    exit
                elseif(this%n(i)<b%n(i))then
                    flag=.false.
                    exit
                end if
            end do
        end if
        ha=max(this%h,b%h)
        if(flag)then
            carry=0
            do i=1,ha
                a=this%n(i)-b%n(i)-carry
                if(a<0)then
                    c%n(i)=(maxnum+1)+a
                    carry=1
                else
                    c%n(i)=a
                    carry=0
                end if
            end do
            c%s=1
        else
            carry=0
            do i=1,ha
                a=b%n(i)-this%n(i)-carry
                if(a<0)then
                    c%n(i)=(maxnum+1)+a
                    carry=1
                else
                    c%n(i)=a
                    carry=0
                end if
            end do
            c%s=-1
        end if
        do i=ha,1,-1
            if(c%n(i)/=0)then
                c%h=i
                return
            end if
        end do
        c%h=1
    end function sub_pos

    function mul_big_big(this,b)result(c)
        class(big_int),intent(in)::this,b
        type(big_int)::c
        integer::ha,hb
        integer::i,j
        if(b==big_int_zero)then
            c=big_int_zero
            return
        end if
        ha=this%h
        hb=b%h
        do i=1,ha
            do j=1,hb
                c%n(i+j-1)=c%n(i+j-1)+this%n(i)*b%n(j)
            end do
        end do
        do i=1,ha+hb-1
            c%n(i+1)=c%n(i+1)+c%n(i)/10
            c%n(i)=mod(c%n(i),10)
        end do
        if(c%n(ha+hb)==0)then
            c%h=ha+hb-1
        else
            c%h=ha+hb
        end if
        c%s=this%s*b%s
    end function mul_big_big

    function mul_big_small(this,x)result(c)
        class(big_int),intent(in)::this
        type(big_int)::b
        integer,intent(in)::x
        type(big_int)::c
        integer::ha,hb
        integer::i,j
        if(x==0)then
            c=big_int_zero
            return
        end if
        b=x
        ha=this%h
        hb=b%h
        do i=1,ha
            do j=1,hb
                c%n(i+j-1)=c%n(i+j-1)+this%n(i)*b%n(j)
            end do
        end do
        do i=1,ha+hb-1
            c%n(i+1)=c%n(i+1)+c%n(i)/10
            c%n(i)=mod(c%n(i),10)
        end do
        if(c%n(ha+hb)==0)then
            c%h=ha+hb-1
        else
            c%h=ha+hb
        end if
        c%s=this%s*b%s
    end function mul_big_small

    function add(this,b)result(c)
        class(big_int),intent(in)::this,b
        type(big_int)::tmpa,tmpb
        type(big_int)::c
        integer::sa,sb
        tmpa=this
        tmpb=b
        sa=tmpa%s
        sb=tmpb%s
        if(sa==1)then
            if(sb==1)then
                !c = a + b
                c=add_pos(tmpa, tmpb)
            else
                !c = a + (-b)
                tmpb%s=1
                c=sub_pos(tmpa, tmpb)
            end if
        else
            if(sb==1)then
                !c = b + (-a)
                tmpa%s=1
                c=sub_pos(tmpb,tmpa)
            else
                !c = (-a) + (-b)
                tmpa%s=1
                tmpb%s=1
                c=add_pos(tmpa, tmpb)
                c%s=-1
            end if
        end if
    end function add

    function sub(this,b)result(c)
        class(big_int),intent(in)::this,b
        type(big_int)::tmpa,tmpb
        type(big_int)::c
        integer::sa,sb
        tmpa=this
        tmpb=b
        sa=tmpa%s
        sb=tmpb%s
        tmpa%s=1
        tmpb%s=1
        if(sa==1)then
            if(sb==1)then
                !c = a - b
                c=sub_pos(tmpa, tmpb)
            else
                !c = a - (-b) = a + b
                c=add_pos(tmpa, tmpb)
            end if
        else
            if(sb==1)then
                !c = (-a) - b =-( a + b)
                c=add_pos(tmpa,tmpb)
                c%s=-1
            else
                !c = (-a) - (-b) = b - a 
                c=add_pos(tmpb, tmpa)
            end if
        end if
    end function sub

    function div(this,b)result(c)
        class(big_int),intent(in)::this,b
        type(big_int)::tmpa,tmpb
        type(big_int)::c
        tmpa=this
        tmpb=b
        tmpa%s=1
        tmpb%s=1
        c=div_pos(tmpa, tmpb)
        c%s=this%s*b%s
    end function div

    function div_pos(this,b)result(c)
        type(big_int),intent(in)::this,b
        type(big_int)::c,tmp
        integer::i,j,ha,hb,hc,hd,ht
        if(b==big_int_zero)then
            error stop "division by zero"
        end if
        if(this<b)then
            c=0
        elseif(this==b)then
            c=1
        else
            ha=this%h
            hb=b%h
            if(this%n(ha)>=b%n(hb))then
                hc=ha-hb+1
                hd=hb
            else
                hc=ha-hb
                hd=hb+1
            end if
            c%h=hc
            tmp%n(1:hd)=this%n(ha-hd+1:ha)
            tmp%h=hd
            do i=1,9
                if(b*i<=tmp.and.b*(i+1)>tmp)then
                    c%n(hc)=i
                    exit
                end if
            end do
            tmp=tmp-b*i
            do i=1,ha-hd
                if(tmp==big_int_zero)then
                    ht=1
                    tmp%h=ht
                else
                    ht=tmp%h
                    tmp%n(2:ht+1)=tmp%n(1:ht)
                    tmp%h=ht+1
                end if
                tmp%n(1)=this%n(ha-hd+1-i)
                do j=0,9
                    if(b*j<=tmp.and.b*(j+1)>tmp)then
                        c%n(hc-i)=j
                        exit
                    end if
                end do
                tmp=tmp-b*j
            end do
        end if
    end function div_pos


    function mod_big_pos(this,b)result(tmp)
        type(big_int),intent(in)::this,b
        type(big_int)::c,tmp
        integer::i,j,ha,hb,hc,hd,ht
        if(b==big_int_zero)then
            error stop "division by zero"
        end if
        if(this<b)then
            tmp=this
        elseif(this==b)then
            tmp=0
        else
            ha=this%h
            hb=b%h
            if(this%n(ha)>=b%n(hb))then
                hc=ha-hb+1
                hd=hb
            else
                hc=ha-hb
                hd=hb+1
            end if
            c%h=hc
            tmp%n(1:hd)=this%n(ha-hd+1:ha)
            tmp%h=hd
            do i=1,9
                if(b*i<=tmp.and.b*(i+1)>tmp)then
                    c%n(hc)=i
                    exit
                end if
            end do
            tmp=tmp-b*i
            do i=1,ha-hd
                if(tmp==big_int_zero)then
                    ht=1
                    tmp%h=ht
                else
                    ht=tmp%h
                    tmp%n(2:ht+1)=tmp%n(1:ht)
                    tmp%h=ht+1
                end if
                tmp%n(1)=this%n(ha-hd+1-i)
                do j=0,8
                    if(b*j<=tmp.and.b*(j+1)>tmp)then
                        c%n(hc-i)=j
                        exit
                    end if
                end do
                tmp=tmp-b*j
            end do
        end if
    end function mod_big_pos

    function mod_big(this,b)result(c)
        !c=this-this/b*b
        class(big_int),intent(in)::this,b
        type(big_int)::tmpa,tmpb
        type(big_int)::c
        tmpa=this
        tmpb=b
        tmpa%s=1
        tmpb%s=1
        c=mod_big_pos(tmpa, tmpb)
        c%s=this%s*b%s
    end function mod_big

    subroutine init_big_int_char(this,b)
        class(big_int),intent(out)::this
        character(len=*),intent(in)::b
        character(len=:),allocatable::s
        integer::m,num,i
        integer(kind=1)::j
        s=trim(adjustl(b))
        m=len_trim(s)
        this%s=1
        if(s(1:1)=="-")then
            m=m-1
            this%s=-1
            s=trim(adjustl(s(2:)))
        elseif(s(1:1)=="+")then
            m=m-1
            s=trim(adjustl(s(2:)))
        end if
        num=m
        do i=0,m-1,1
            read(s(i+1:i+1),*)j
            this%n(num-i)=j
        end do
        this%h=num
    end subroutine init_big_int_char

    subroutine init_big_int_int(this,x)
        class(big_int),intent(out)::this
        integer,intent(in)::x
        character(len=13)::str
        str=""
        write(str,"(g0)")x
        this=str
    end subroutine init_big_int_int

    subroutine write_big_int(this)
        class(big_int)::this
        integer::m
        integer::i
        m=this%h
        write(*,"(I2)",advance="no")this%n(m)*this%s
        do i=m-1,1,-1
            write(*,"(I1)",advance="no")this%n(i)
        end do
        write(*,*)
    end subroutine write_big_int


end module

program main
    use big_int_mod
    implicit none
    type(big_int)::a,b,c
    integer(kind=8)::a1,b1,c1
    a="12345678901234567"
    a1=12345678901234567_8
    call a%write_big_int()
    b="987654321"
    b1=987654321
    call b%write_big_int()
    write(*,*)"c=a+b",a1+b1
    c=a+b
    call c%write_big_int()
    write(*,*)"c=a-b",a1-b1
    c=a-b
    call c%write_big_int()
    write(*,*)"c=a*b"
    c=a*b
    call c%write_big_int()
    write(*,*)"c=a/b",a1/b1
    c=a/b
    call c%write_big_int()
    write(*,*)"c=mod(a,b)",mod(a1,b1)
    c=mod(a,b)
    call c%write_big_int()
    a="123456789012345678901234567890123456789"
    write(*,*)"a="
    call a%write_big_int()
    write(*,*)"b="
    b="98765432109876543210"
    call b%write_big_int()
    write(*,*)"c=a+b"
    c=a+b
    call c%write_big_int()
    write(*,*)"c=a-b"
    c=a-b
    call c%write_big_int()
    write(*,*)"c=a*b"
    c=a*b
    call c%write_big_int()
    write(*,*)"c=a/b"
    c=a/b
    call c%write_big_int()
    write(*,*)"c=mod(a,b)"
    call write_big_int(a-c*b)
    c=mod(a,b)
    call c%write_big_int()
end program main
